How synaptic function controls critical transitions in spiking neuron networks: insight from a Kuramoto model reduction

The dynamics of synaptic interactions within spiking neuron networks play a fundamental role in shaping emergent collective behavior. This paper studies a finite-size network of quadratic integrate-and-fire neurons interconnected via a general synaptic function that accounts for synaptic dynamics and time delays. Through asymptotic analysis, we transform this integrate-and-fire network into the Kuramoto-Sakaguchi model, whose parameters are explicitly expressed via synaptic function characteristics. This reduction yields analytical conditions on synaptic activation rates and time delays determining whether the synaptic coupling is attractive or repulsive. Our analysis reveals alternating stability regions for synchronous and partially synchronous firing, dependent on slow synaptic activation and time delay. We also demonstrate that the reduced microscopic model predicts the emergence of synchronization, weakly stable cyclops states, and non-stationary regimes remarkably well in the original integrate-and-fire network and its theta neuron counterpart. Our reduction approach promises to open the door to rigorous analysis of rhythmogenesis in networks with synaptic adaptation and plasticity.


Introduction
Cooperative rhythms play a pivotal role in brain functioning.Fully or partially synchronized oscillations, observed across various frequency bands, underlie fundamental processes such as perception, cognition, and motor control Churchland and Sejnowski, 1992;Mizuseki and Buzsaki, 2014;Kopell et al., 2000.Extensive research has focused on the emergence of cooperative rhythms in networks of spiking and bursting neurons, encompassing synchronization Kopell et al., 2000;Brunel, 2000;Börgers and Kopell, 2003;Somers and Kopell, 1993;Izhikevich, 2007;Belykh et al., 2005;Ermentrout and Terman, 2010, partial and cluster synchronization Achuthan and Canavier 2009;Shilnikov et al., 2008;Belykh and Hasler, 2011;Schöll, 2016;Berner et al., 2021a, neural bumps Laing andChow, 2001;Gutkin et al., 2001, andchimera states Olmi et al., 2011;Omelchenko et al., 2013.Networks of spiking neurons with fast synaptic connections are often modeled via pulsatile on-off coupling, which sharply activates upon the arrival of a spike from a pre-synaptic cell.Such interactions are conveniently represented by networks of quadratic integrate-and-fire (QIF) models particularly suitable for large-scale simulations and analysis of cooperative dynamics Gerstner and Kistler (2002).The macroscopic dynamics of QIF networks have received extensive attention through the reduction to low-dimensional model descriptions, especially in the thermodynamic limit of infinitedimensional networks Montbrió et al., 2015;Pazó and Montbrió, 2016;Devalle et al., 2017;Esnaola-Acebes et al., 2017;Devalle et al., 2018;Schmidt et al., 2018;Pietras et al., 2019;Montbrió and Pazó 2020;Lin et al., 2020;Gast et al., 2020;Taher et al., 2020;Byrne et al., 2022;Clusella et al., 2022;Clusella and Montbrió 2024;Ratas and Pyragas, 2018;Pyragas andPyragas 2022, 2023;Coombes 2023;Ferrara et al., 2023. Notably, Montbrió et al. (2015) derived exact macroscopic equations for QIF networks, uncovering an effective coupling between firing rate and mean membrane potential governing network dynamics.Pietras et al. (2023) offered an analytical description of QIF network macroscopic dynamics, extending beyond the Ott-Antosen ansatz Ott and Antonsen (2008) and exploring various fast synaptic pulse profile choices.The impact of synaptic time delay on the collective dynamics of integrate-and-fire networks with sharply activated synaptic coupling, modeled by the Dirac delta function, was also extensively explored (Ernst et al., 1995;Devalle et al., 2018;Ratas and Pyragas, 2018;Pyragas andPyragas 2022, 2023).In particular, Devalle et al. (2018) reduced a QIF model with synaptic delay to a set of firing rate equations to analyze the existence and stability of partially synchronous states.Ratas and Pyragas (2018) employed a Lorenzian ansatz to characterize macroscopic oscillations of a QIF network with heterogeneous time-delayed delta function synapses.However, there is a lack of analytical studies on the role of slower synaptic activation, potentially in the presence of time delays, in controlling critical phase transitions in QIF networks.Nevertheless, since the seminal paper by Van Vreeswijk et al. (1994), it has been recognized that slow inhibitory and excitatory synapses can reverse their roles, with slow inhibition favoring synchronization Golomb and Rinzel 1993;Terman et al., 1998;Elson et al., 2002.While predicting the exact rates of synaptic activation inducing such critical transitions in conductance-based spiking models may be challenging, analytically tractable QIF networks offer promising avenues for such exploration.
Toward this goal, this paper investigates a finite-size network of QIF neurons globally connected via a general kernel function that governs both synaptic activation and synaptic time delay.We analytically illustrate how the shape of the kernel function impacts neuron interaction, significantly altering the microscopic and macroscopic behavior of QIF networks representing Type I neuron populations.This is achieved by reducing QIF networks and their phase analog, theta neuron networks, to the Kuramoto-Sakaguchi (KS) model.Here, oscillator frequencies, coupling strength, and the Sakaguchi phase lag parameter are determined by the pulse profile's first and second terms in the Fourier expansion.We conduct this reduction under the weak coupling assumption, utilizing the intermediate step of representing the QIF network as a generalized Winfree model, subsequently reduced to the KS model.
In our recent study Munyayev et al. (2023), we elucidated the qualitative connection between the dynamics of QIF networks incorporating synaptic dynamics and neuronal refractoriness, and the second-order Kuramoto model with high-order mode coupling.Here, we use multi-scale analysis to derive exact relationships between the QIF network with an arbitrary synaptic activation function and the KS model.Specifically, we establish explicit conditions on the parameters of the general kernel function that lead to critical transitions, determining whether the coupling is attractive or repulsive.Consequently, these conditions dictate the emergence of stable synchronization or nonstationary generalized splay states Berner et al. (2021b) and cyclops states Munyayev et al. (2023).Our analysis reveals alternating stability regions for network synchronization, dependent on both the (slow) synaptic activation and time delay.With some important caveats, this finding can be interpreted as an analogous stability criterion for synchronization in timedelayed phase oscillator networks Earl and Strogatz (2003).
Our approach serves as a connecting link between two alternative methodologies for describing macroscopic dynamics: QIF networks and theta neurons, and Winfree-type models Pazó and Montbrió, 2014;Gallego et al., 2017;Montbrió and Pazó, 2018;Pazó et al., 2019;Pazó and Gallego 2020;Manoranjani et al., 2021;Bick et al., 2020.Our KS model reduction of the generalized Winfree model with a general synaptic activation function can be seen as an extension of the work Montbrió and Pazó (2018), where a twopopulation Kuramoto model was derived from a network of Winfree oscillators featuring a feedback loop between fast excitation and slow inhibition.
The structure of this paper is outlined as follows.Section 2 presents the QIF network model, its theta neuron equivalent, and the general synaptic activation function.Section 3 details transforming the theta neuron model into the generalized Winfree model.We expand the pulse profile as a Fourier series and further simplify the model to the KS model using weak coupling-enabled averaging techniques.Section 4 focuses on a specific example of synaptic activation, presenting a class of kernel functions.We establish exact conditions determining whether the synaptic coupling is attractive, promoting synchronization, or repulsive, favoring splay and cyclops states.Section 5 offers numerical validation of the derived conditions and presents a comparison between the dynamics of the QIF network, the theta neuron model, and the reduced KS model.We demonstrate that the KS model accurately predicts firing rates and times, capturing the emergence of synchronization, weakly stable cyclops states, and non-stationary regimes.Section 6 contains concluding remarks and discussions.

The general QIF network and its theta neuron representation
Physiologically, excitable neurons are commonly categorized into two types.We focus here on Type I neurons, a group encompassing cortical excitatory pyramidal neurons.When subjected to a sufficiently large input stimulus, these neurons exhibit action potentials at an arbitrarily low rate, signaling the disappearance of a resting state through a saddle-node bifurcation.The canonical model used to describe Type I neurons is the QIF neuron model, which characterizes neurons' dynamics near the spiking threshold Izhikevich (2007).
This study investigates a globally coupled network of N QIF neurons interacting through chemical synapses.Each neuron's microscopic state is characterized by its individual membrane potential v n , governed by the following ordinary differential equation Izhikevich (2007): Here, η n , n 1, 2, . . ., N represents external constant currents applied to neurons, ϰ is a common synaptic weight controlling the total strength of synaptic inputs, and S(t) is a time-varying input drive.When the membrane potential v n of the nth neuron reaches the threshold value v th , the neuron generates a spike, and its voltage resets to v r .In the absence of the input drive (S(t) 0), the intrinsic applied current η n 0 places the corresponding nth neuron at a saddle-node bifurcation, marking the onset of periodic firing.Thus, if η n < 0, the neuron is in the excitable regime, while if η n > 0, it is in the oscillatory regime.
The last term on the right-hand side of (Eq. 1) represents synaptic interactions characterized by the coupling strength ϰ of the global synaptic drive.In general, the mean population synaptic activity S(t) can be expressed by the following recurrent input equation: This equation accounts for relaxation processes and describes a specific type of neuron activation and its sensitivity to stimuli from other cells, including signal duration and post-spike latency.Here, t m n′ denotes the time of the mth spike of the n′th neuron, δ(t) represents the Dirac delta function, and G(t/τ) is the normalized synaptic activation caused by a single presynaptic spike with a time scale τ.Notably, the integral transformation with the kernel G(t/τ) acts as a low-pass filter.
The QIF-neuron model (Eq. 1) describes the membrane potential v n and operates as a hybrid dynamical system, incorporating instantaneous resets to a base value v r upon spike emission.While this formulation provides a direct physical interpretation, discontinuities can pose challenges for certain applications.Fortunately, a smooth change of coordinates exists, transforming the QIF-neuron dynamics into a space where the membrane potential v n is represented by a phase variable θ n on the unit circle.This representation captures nonlinear spikegenerating mechanisms of Type I neurons, ensuring smooth solutions within a compact domain.In the limit v th −v r → ∞, the transformation v n (t) tan(θ n (t)/2) Ermentrout and Kopell (1986) converts the membrane potential description of the QIFneuron model (Eq. 1) into a canonical theta neuron model of a population of Type I neurons coupled by an excitatory or inhibitory synaptic drive.In this case, each neuron's dynamics is governed by the following equation: where n 1, . . ., N represents the index of the nth neuron, and its state is characterized by the phase angle θ n .We assume that the constant excitability parameters η n , akin to fixed input currents, Synaptic dynamics S(t), defined by (Eq.4), induced by presynaptic spikes p(t), taking the form of (A,B) P ∞ (t) as defined by (Eq. 6) and (C,D) P ] (t) as defined by (Eq.5) with ] 40 and linearly increasing phase θ 1 (t) t (N 1).The parameter τ 0.18 is used for all cases.
have slightly different values across the network elements.Moreover, we consider η n > 0, placing each neuron in an oscillatory regime indicative of periodic spiking.A neuron is deemed to spike or produce an action potential when θ n crosses π while increasing.Consequently, in addition to a common external input η n , each cell receives stimulation from other cells.Thus, the neurons are recurrently coupled via synaptic current S(t).
The last term on the right-hand side of (Eq. 3) accounts for chemical interactions among neurons.The coupling strength ϰ is assumed to be uniform across all neurons.We model the synaptic activity S(t) acting on a neuron with the following expression: where the function determines the shape of the pulsatile chemical synapse.The positive integer parameter ] controls the sharpness of P ] (θ n ), with higher values yielding sharper peaks.Note that as ] → ∞, the smooth profile P ] (θ n ) converges to P ∞ (θ n ), effectively representing δ-pulses so that where C T (•) denotes the Dirac comb with period T. In this limit and under the assumption v th −v r → ∞, the theta neuron model (Eqs 3, 4) fully aligns with the QIF-neuron model (Eqs 1, 2).It is noteworthy that these models exhibit an unconventional characteristic: in contrast to conductance-based models, inhibitory δ-pulse coupling (ϰ < 0) promotes synchronization, thus being considered attractive, while excitatory δ-pulse coupling (ϰ > 0) is repulsive Izhikevich ( 2007).
Note that the limiting δ-pulse case with P ∞ (θ n ) offers a convenient framework in which function G(t) controls synaptic activation and deactivation after the instantaneous occurrence of a presynaptic action potential.The use of P ∞ (θ n ) simplifies the analytical expressions in the KS model to be more explicit and manageable (see next section).However, the general Eqs 4, 5 for synaptic input S(t) provide a more comprehensive and accurate description that captures the nuances of synaptic transmission processes.The release of neurotransmitters from the presynaptic neuron that induces a chemically activated synaptic current is non-instantaneous.Thus, the general function P ] (θ n ) more adequately describes this gradual neurotransmitter release that takes place as the presynaptic neuron state approaches a spike activation threshold.
In this work, we primarily focus on the pulse shape defined by (Eq.5), originally proposed in Ariaratnam and Strogatz (2001) and widely adopted in recent studies of pulse-coupled phase oscillators O'Keeffe and Strogatz 2016; Pazó and Montbrió, 2014;Gallego et al., 2017;Pazó et al., 2019;Bick et al., 2020 and populations of theta neurons Luke et al., 2013;So et al., 2014;Laing 2014, Laing 2015;Montbrió et al., 2015;Pazó and Montbrió, 2016;Chandra et al., 2017;Goel and Ermentrout 2002;Bick et al., 2020;Pietras et al., 2023.However, our approach is directly applicable to alternative pulse shapes satisfying common properties such as unimodality, normalization, symmetry, and localization around θ n π and considered in previous studies Gallego et al., 2017;Pietras et al., 2023.In the following, we explore how the shape of the kernel function G(t/τ), governing synaptic activation, impacts the network's collective behavior.To tackle this analytically, we will demonstrate that the model (Eqs 3, 4) consideration can be reduced to a more analytically-tractable KS model.The process involves several key steps: firstly, leveraging the assumption η n > 0, we introduce an alternative phase representation for (Eqs 3, 4).Subsequently, we will derive the KS model of phase oscillators by employing multiple time scale analyses in weak chemical synaptic coupling scenarios.

FIGURE 2
Characteristics of the synaptic dynamic profile for S(t) (Eq.4) as a function of q. (A) Peak latency, t m , (B) maximum value, S m , and (C) full width at half maximum (FWHM), induced by spikes P ] (t) with τ 0.1 and different ] (cyan markers -] 10, red markers -] 100, orange markers -] 1000).

Deriving the KS model from the thetaneuron model: an asymptotic analysis
We begin by assuming that each neuron operates within an oscillatory regime in the absence of interaction, i.e., η n > 0. Hence, we can use a dynamical variable transformation: where Ω is an unknown parameter to be determined subsequently.It is notable that Ω is close to 2 〈η n 〉 , where 〈η n 〉 denotes the mean value of the distributed external constant currents η n .However, a more precise determination of Ω is feasible for finite-size networks, as will be shown later in this section.
The transformation (Eq.7) transitions the model (Eqs 3, 4) to an alternative phase representation: where This representation remains consistent with the original description; notably, φ n ∈ [−π, π], and the occurrence of spikes for the nth neuron is Regions of attractive (blue) and repulsive (red) coupling in the theta neuron model (Eqs 3, 4), corresponding to the KS model regions defined by (Eq.29).The colors represent the coupling strength K sign as a function of the synaptic time constant τ and the common external input θ for a fixed q. (A) q 2, (B) q 4, and (C) q 12.Other parameters: ϰ −0.2π, ] 20, η 1 η 2 / η N η.The yellow points A, B, C, and D indicate the parameter values used for numerical simulations of Figures 6-9.
Coupling strength K (A) and Sakaguchi parameter α (B) as functions of synaptic adaptation/delay (parameter q) and the finite pulse width (parameter ]).The dependencies are calculated analytically via (Eq.21) and (Eq.27) for τ 0.15, ϰ −0.2π, η 2. The coupling strength K increases as the pulse width decreases (via increasing ]) and decreases as synaptic adaptation slows down and experiences larger time delays (via increasing q).The Sakaguchi parameter α is highly sensitive to q and only weakly dependent on ].
still defined by φ n crossing π.However, given that in the model (Eqs 3, 4), all cells receive constant external inputs η n > 0, and all units operate within the oscillatory regime, each phase φ n uniformly rotates in the absence of network interaction.Hence, the representation (Eqs 8, 9) proves more convenient for subsequent analysis.Thus, we derive the governing equation for the introduced dynamical variables φ n that may be viewed as the generalized Winfree model, which, in turn, can be reduced to the KS model.Further elaboration on this approach's derivation and technical intricacies are presented below.
For further analysis, it is convenient to expand the symmetric pulse Q ] (φ n ) into a Fourier series with respect to φ n .This expansion takes the following form: The coefficients Q ]ℓ in this series can be expressed analytically as follows: where C 2m 2ℓ represents combinations, Γ(•) denotes the gamma function, and 2 F 1 (•, •; •; •) is the hypergeometric function.Specifically, for the first two Fourier series coefficients Q ∞0 and Q ∞1 of the symmetric pulse Q ∞ (φ n ) in (Eq.10), we obtain: Noteworthy, in the limit ] → ∞, i.e., for Q ∞ (φ n ), all coefficients in the Fourier series (Eq.11) converge to the same value Q ∞ℓ (−1) ℓ |Ω|/2π.We proceed by assuming that the synaptic coupling is weak, allowing us to express it as 2ϰ/Ω εκ, where ε ≪ 1 is a small parameter.Similarly, we assume that the deviations of external inputs η n from the value Ω 2 /4 are small, i.e., η n − Ω 2 /4 εΩσ n /2.
These assumptions enable a multiple-time scale analysis.To facilitate this analysis, we introduce a separation of time scales: and represent each phase variable, φ n (t), as an asymptotic series with respect to the small parameter ε: Substituting the series (Eq.15) with times (Eq.14) and (Eq.11) into (Eqs 8, 9), and considering the zeroth order in ε, we obtain for φ (0)  n (t 0 , t 1 , t 2 , . ..): and, taking into account (Eq.16), for S(t 0 , t 1 , t 2 , . ..) we arrive at where each corresponding complex coefficient G ℓ is determined as follows: In (Eq.16), the first term Ωt 0 describes the fast, free-running rotation of period 2π/Ω, while the slow phase drifts induced by synaptic interaction are characterized by the set of slow variables ϕ n (t 1 , t 2 , . ..).Hence, ϕ n (t 1 , t 2 , . ..) can be considered constant over time scales comparable to the period of the corresponding fast rotation.Consequently, the standard averaging method can be applied to derive the KS model corresponding to (Eqs 3, 4).Notably, in this case, the Sakaguchi phase shift emerges naturally due to the complex nature of the coefficient G ℓ determined in (Eq.18).
In accordance with the averaging procedure after substituting expression (Eq.17) into (Eq.8), the next step of our asymptotic approach involves considering all terms that are o(ε) and, in the first order in ε, we obtain a set of equations for φ (1)  n (t 0 , t 1 , . ..).To eliminate the secular terms that grow without bounds as t 0 → ∞, we impose the conditions This yields a solution for φ (1) n (t 0 , t 1 , . ..) without the secular terms.Note that (Eq.19) measures the rate of change of ϕ n (t 1 , t 2 , . ..) with respect to the slow time scale t 1 .Finally, taking into account the relation dϕ n /dt ≈ ε∂ϕ n /∂t 1 , we find that the dynamics of the slow phases ϕ n (t) is approximately described by the KS model: where To determine the unknown parameter Ω, we set the mean value of the intrinsic frequencies ω n in the KS model (Eq.20) to zero.Hence, the value of Ω can be found by solving the nonlinear algebraic equation: This choice of the optimal value of parameter Ω yields a better quantitative match between the numerical simulation results of the theta neuron model (Eqs 3, 4) and the KS model (Eq.20), compared to the conventional choice of 2 〈η n 〉 .In the limiting case ] → + ∞, where the shape of the pulsatile chemical synapse P ∞ (θ n ) is determined by Eq. 6, the first two Fourier series coefficients Large-size networks.Dynamical equivalence between the theta neuron (Eqs 3, 4) and KS models (Eq.20), demonstrated via the onset of full synchronization (A,B) and non-stationary generalized splay state with an oscillating |R 1 |≈ 0 (C,D) for N 100 (A,C) and N 1, 000 (B,D).Notations are as in Figures 5, 6. Parameters: q 2, ϰ −0.2π, τ 0.5, ] 20, η 2.0 (A,B); q 2, ϰ −0.2π, τ 0.8, and ] 20 (C,D).
Eqs 8, 9 converge to Q ∞0 −Q ∞1 |Ω|/2π and, hence, Eq. 22 for Ω can be solved analytically.This yields the simpler expressions for the parameters ω n , K and α of the KS model (Eq.20) with coefficients (Eqs 21a, 21b, 21c): where the complex coefficient G 1 is determined as follows: Note that oscillator frequencies ω n , coupling strength K, and the Sakaguchi phase lag parameter α in the KS model (Eq.20) are explicitly defined through by the pulse profile's first and second Fourier series terms Q ]0 and Q ]1 .In particular, the expressions (Eq.23) clearly indicate that for the δ-pulses, the sign of the coupling strength K in the KS model is solely determined by and is opposite to the sign of the coupling κ in the original QIF model.The following section will show that this property carries over to the general finite-width pulses defined by (Eq.5).For a specific class of the activation function G(t/τ), we will also demonstrate that increasing the pulse width decreases the coupling strength K and practically does not affect the Sakachichi phase lag parameter α.The direct dependence on the properties of synaptic activation enables a straightforward assessment of the role of synaptic interactions in critical phase transitions and dynamics.Specifically, the coupling strength K and Sakachichi parameter α directly reflect the impact of synaptic activation on critical phase transitions and the synchronization behavior of the neuron population.In particular, determining whether the coupling in the theta neuron model (Eqs 3, 4) is attractive or repulsive presents a challenge due to the complexity of the system, whereas the KS model (Eq.20) makes this process straightforward.In the subsequent section, we delve into this process, focusing on a particular synaptic activation profile.Onset of full synchronization in the QIF (Eq. 1) and KS models (Eq.20).(A) Firing rate and (B) firing times of QIF neurons (cyan curves and round markers) and oscillators of the KS model (black curves and cross markers) with the firing times recorded at θ n (t f ) π.Each row in (B) represents the firing times of a neuron/oscillator.Inset (C) zooms-in on the firing time pattern from (B).Initial conditions are as in Figure 5. Parameters: N 21, q 4, τ 0.15, ϰ −0.2π, v th −v r 10 5 , η 2, δη 6 × 10 −3 , ] 10 5 correspond to point C in Figure 3.The firing rate was calculated within a sliding time window of δt 5 × 10 −2 .

The role of synaptic profile: a combined effect of activation, deactivation, and time delays
To demonstrate the important effects arising from the specific selection of the shape of a "low-pass filter" in synaptic activation, we examine the following class of kernel functions: where H(t/τ) is the Heaviside step function and q represents an integer parameter pivotal to the network dynamics, as will become evident subsequently.The chosen kernel corresponds to the Green's function for a non-homogeneous linear differential equation of order (q + 1).In this case, the equation for the mean synaptic activity S(t) is generated by the (q + 1)th order linear differential operator L (τ d/dt + 1) q+1 and contains an external signal that represents an average profile of spike pulses: ), The source term on the right-hand side of (Eq.26) is presented in three interchangeable forms corresponding to the QIF model (Eqs 1, 2), theta neuron (Eqs 3, 4), and their averaged representation through the KS model (Eq. 20).This term can be interpreted as the population firing rate, which induces a post-synaptic current in response to the arrival of spikes.In the limit τ → 0, one might assume that the interaction between neurons becomes instantaneous.However, if the characteristic time scale τ of a post-synaptic response is not negligibly small, it becomes imperative to take into account the individual adaptation dynamics of the synaptic variable S(t), encompassing its activation, deactivation, and time delay.
In the case q 0, the mean population synaptic activity S(t) is governed by the standard relaxation rule.Here, when the n′th neuron fires at time t m n′ , generating the mth spike in the form of the Dirac delta function, the mean activity S(t) instantaneously changes and subsequently decays exponentially in the absence of further firings.The parameter τ acts as a synaptic time constant.
For q 1, when the n′th neuron fires at time t m n′ and the mth Dirac delta pulse is generated, the variable S(t) is augmented by the function G((t − t m n′ )/τ)/N defined by (Eq.25) with q 1, coinciding with the so-called alpha-function pulse Mohanty and Politi 2006;Zillmer et al., 2007;Bolotov et al., 2016;Chen et al., 2017.For this alpha-pulse created by a spiking neuron, τ determines both the signal's width and the time at which it attains its maximum value (Figure 1).
We extend the argument to arbitrary q and τ that make the model of synaptic adaptation (Eq.26) encompass a broad spectrum of realistic biophysical scenarios, ranging from fast non-delayed to slow delayed activation.These scenarios include neurotransmitter release in the synaptic cleft and the opening/closing of postsynaptic ion channels, characterized by distinct time scales such as latency Diagram similar to Figure 8 showing a nearly perfect match for asynchronous firing rate (A) and firing times (B) in the QIF network (cyan curves and round markers) and the KS model (black curves and cross markers).Parameters: N 21, q 4, τ 0.41, ϰ −0.2π, v th −v r 10 5 , η 2, δη 6 × 10 −3 , ] 10 5 correspond to point D in Figure 3.Other notations and settings are as in Figure 8.
Frontiers in Network Physiology frontiersin.org10 (time delay), rise, and decay times, as reflected in the kernel function (Eq.25).
Our choice of the kernel function G(t/τ) for synaptic activation aligns with the gamma distribution Zwillinger (1992), a continuous probability distribution characterized by two parameters: τ > 0, the scale parameter, and q > 0, the shape parameter.This distribution reaches its maximum value at t max qτ, with mean value t mean (q + 1)τ, variance σ 2 t (q + 1)τ 2 , and skewness γ t 2/ q + 1 .The skewness, reflecting the symmetry of the distribution about its mean, is maximal for the exponential case q 0 and diminishes for larger values of q, indicating increased symmetry.While both q and τ influence synaptic dynamics, q has a more significant impact on synaptic time delay than τ, whereas τ predominantly shapes the synaptic activation profile.Figures 1, 2 demonstrate how the parameters τ and q determine the time profile of the post-synaptic response and its characteristics.
Towards our objective of deriving explicit conditions for the attractiveness or repulsiveness of synaptic coupling governed by (Eq.24) with the kernel function (Eq.25), we calculate the complex coefficient G 1 and its modulus and argument from (Eq. 24) as follows: where Ω is determined from (Eq. 22).While the hypergeometric function Q ]1 , a factor defining the coupling strength K in the KS model (Eq.20) with coefficients (Eq.27), cannot be expressed via elementary functions (except for the liming δ-pulse case with ] → ∞), it is evident that Q ]1 ≤ 0 for ] ≥ 0. Consequently, the coupling strength K 2ϰ|G 1 |Q ]1 /Ω in the KS model (Eq.20) for Ω > 0 and the coupling strength κ in the QIF model (Eq. 1) have opposite signs.Therefore, for ϰ < 0, a positive coupling strength K corresponds to attractive coupling, provided that the Sakaguchi phase lag parameter α < π/2.According to (Eq.21), cos α > 0 if Solving the inequality (Eq.28), we obtain the following q/2 + 1 intervals of parameters that correspond to the attractive coupling in the KS model (Eq.20) and therefore in the QIF (Eqs 1, 2) and thetaneuron models (Eqs 3, 4) with κ < 0: Firing rate (A) and firing times (B) of a three-cluster cyclops state in the QIF network (cyan curves and round markers) and the KS model (black curves and cross markers).(C) Snapshots of the cyclops state phase distributions φ n in the KS model at two time instants.The oscillators' coloring represents their phase.The cyclops states is formed by a solitary oscillator (red) and two coherent clusters, each composed of 10 oscillators (orange and blue).The initial phases are chosen near a cyclops state.Parameters: N 21, q 2, τ 0.8, ϰ −0.2π, v th −v r 10 5 , η 2, δη 0, ] 10 5 .Other notations and settings are as in Figure 8.
Figure 3 displays the regions, as defined by (Eq.29), where the coupling in the KS model is attractive (blue) or repulsive (red).These regions exhibit an alternating pattern as functions of the synaptic time constant τ and the common external input θ for a fixed q.Increasing the parameter q, which primarily controls the time delay, makes the synaptic coupling more sensitive and Firing rate (A) and firing times (B) in the QIF network demonstrating the transition to full synchronization starting from random initial conditions.Parameters N 1, 000, q 4, τ 0.15, ϰ −0.2π, v th −v r 100, η 2, and δη 6 × 10 −3 correspond to point C in Figure 3.The simulations use the algorithm from Pazó and Montbrió (2016), which accounts for the neurons' refractory time.The mean synaptic activation was calculated within a sliding time window of δt 5 × 10 −2 .

FIGURE 12
Firing rate (A) and firing times (B) in the QIF network accompanying the transition from synchronous to asynchronous dynamics.Parameters N 1, 000, q 4, τ 0.41, ϰ −0.2π, v th −v r 100, η 2, and δη 6 × 10 −3 correspond to point D in Figure 3.Other notations and settings are as in Figure 11.
results in more alternating zones.This alternating pattern of attractive and repulsive coupling resembles the stability criterion for synchronization in time-delayed phase oscillator networks, where the time delay controls the sign of the derivative of the periodic coupling function Earl and Strogatz (2003).
Figure 4 provides further insight into how the synchronization properties of the KS model, critically controlled by the coupling strength K and Sakaguchi parameter α, depend on the original parameters of the QIF model.Notably, both K and α exhibit a weak dependence on ], suggesting that the pulse width has minimal impact on synchronization dynamics, except for the case of biologically irrelevant wide pulses where ] approaches 0, making the coupling strength K significantly weaker.In contrast, both K and α are highly sensitive to variations in q, with the latter producing the alternating regions of attractive and repulsive coupling, as illustrated in Figure 3.
In the following, we offer additional evidence supporting the predictive power of the derived KS model.We show numerically that it effectively captures the emergence of robust dynamical regimes like synchronization and more intricate partially synchronized dynamics such as weakly stable cyclops states and non-stationary generalized splay states in both the QIF and theta neuron models.

Dynamical equivalence of the models: numerical validation
We conduct numerical computations using a widely accepted fifth-order Runge-Kutta method with a fixed time step of 0.01, providing additional validation for our analytical findings and predictions.
To characterize the dynamical regimes, we utilize both microscopic measures (pairwise phase differences and firing times) and macroscopic indicators such as the firstand second-order complex Kuramoto parameters Daido 1992;Skardal et al., 2011: where r ℓ and ψ ℓ , l 1, 2 define the magnitude and the phase of the ℓth moment Kuramoto order parameter R ℓ (t), respectively.The first-order scalar parameter r 1 |R 1 | characterizes the degree of phase synchrony with r 1 1 corresponding to full phase synchrony.
The second-order scalar parameter r 2 |R 2 | determines the degree of cluster synchrony, where r 2 0 corresponds to generalized splay states Berner et al. (2021b) and their particular case of cyclops states Munyayev et al. (2023).We will also use the firing rate, which measures the average rate at which neurons emit spikes, as another important macroscopic observable to characterize dynamical regimes.We will use the following formula Pietras et al. (2023) in simulations later in the section: To identify the time steps corresponding to neuron spike events in both the QIF-neuron model and the theta neuron model, we monitored the sign changes of v n (t) − v th and (θ n (t) mod 2π) − π, along with their time derivative signs.Subsequently, we determined the spike moment t m n using linear interpolation within each time step.
Figures 5, 6 illustrate the perfect correspondence between the emergence of full synchronization and non-stationary generalized splay states in the theta neuron model (Eqs 3, 4) and the KS model (Eq.20) within the range of attractive coupling (point A in Figure 3) and repulsive coupling (point B in Figure 3), respectively.In the case of full synchronization (Figure 5), the first-order and second-order scalar parameters (Eq.30), |R 1 | and |R 2 |, converge to unity but cannot reach 1 due to intrinsic parameter mismatch in η n .Likewise, the first-order scalar parameter, |R 1 |, associated with the nonstationary generalized splay state oscillates closely around 0 (Figure 6). Figure 7 demonstrates that the KS model maintains its excellent predictive power for the emergence of full synchronization and non-stationary generalized splay states in large networks of 100 and 1,000 theta neurons.
Figures 8-10 illustrate the remarkable agreement between cooperative dynamics in the QIF and KS models.Specifically, Figure 8 depicts the onset of full synchronization, as evidenced by synchronized firing rates and times determined via (Eq.31).The slight discrepancy in the firing times between the QIF and KS models may stem from various sources, such as accumulated numerical errors and the approximate calculation of the frequency parameter Ω derived from (Eq. 22) for selecting the KS model parameters.Figure 9 provides evidence for the capability of the KS model to perfectly predict even nonstationary, asynchronous firing in the QIF model.Figure 10 illustrates the predictive power of the derived KS model in discerning stable complex cluster patterns like cyclops states in the QIF model.Introduced in Munyayev et al. (2023), cyclops states are formed by two distinct, coherent clusters, and a solitary oscillator reminiscent of the Cyclops' eye.While detecting stable cyclops states can be challenging in the QIF model, the KS model provides a more convenient and constructive approach.
In the numerical calculations of relatively small-size networks of QIF neurons presented in Figures 8-10, we employed the fourth-order Runge-Kutta method using the procedure for identifying the spike events described above.To ensure better consistency between the QIF-neuron model and the theta neuron model, we used sufficiently large values of v th and v r , specifically v th −v r 10 5 .We did not account for the time interval required for the membrane potential to reach ± ∞.However, this approach becomes computationally expensive for simulating large networks of QIF neurons.Therefore, in the numerical calculations presented in Figures 11, 12, we employed the algorithm described in Pazó and Montbrió (2016).This algorithm, based on the Euler method, accounts for the time it takes for the dynamic variable associated with the membrane potential to pass through the singularity ± ∞ after exceeding v th and reach the final value v r .A detailed description of the algorithm and its parameters is available in the Supplemental Material for Pazó and Montbrió (2016).Figures 11, 12 confirm that the main analytical predictions for the dynamics of the QIF model remain valid for large network sizes and regardless of the calculation scheme used.Specifically, Figure 11 validates the prediction of coupling attractiveness at point C in Figure 3 and demonstrates the onset of full synchronization in the QIF network of 1,000 neurons.Similarly, Figure 12 confirms the coupling repulsiveness at point D in Figure 3 and illustrates the emergence of asynchronous dynamics.Notably, the asynchronous firing rate in the 21-node QIF network shown in Figure 9, corresponding to point D, may resemble weak coherence patterns.In contrast, its 1,000-node counterpart in Figure 12, with parameters also corresponding to point D, exhibits nearly perfect asynchronous dynamics.Our numerous additional simulations of 1,000-node QIF networks for a large set of parameters, including those near the boundary of the alternating regions of attractive and repulsive coupling in Figure 3, further validated the consistency of the predictions for cooperative dynamics and the critical transitions in the QIF networks [not shown].

Conclusions
Understanding the influence of synaptic dynamics, including activation rates, deactivation processes, and latency, on collective dynamics in neuronal networks is of significant importance.Considerable advancements have been made in analyzing the role of fast or time-delayed synapses in integrate-and-fire neuron networks.However, there remains a scarcity of analytical studies exploring the influence of slower synaptic dynamics, potentially in the presence of time delays, on controlling critical phase transitions in neuronal networks.
In this paper, we have made substantial contributions to advancing analytical methods in this field.We studied a finitesize network of QIF neurons globally interconnected via a generalized kernel function governing both synaptic activation and time delay.Our analytical exploration demonstrated how the shape of the kernel function profoundly affects neuron interaction, thereby significantly modifying the microscopic and macroscopic behavior of QIF networks.To achieve this, we reduced the QIF and theta neuron network models to the Kuramoto-Sakaguchi model.In this model, oscillator frequencies, coupling strength, and the Sakaguchi phase lag parameter are determined by the Fourier terms of the pulse profile series expansion.
We established exact conditions determining whether synaptic coupling is attractive, fostering synchronization, or repulsive, promoting splay and cyclops states.Furthermore, we demonstrated a remarkable correspondence between the dynamics of the derived KS model and the original QIF and theta neuron models.Specifically, the KS model accurately predicted firing rates and times, capturing the emergence of synchronization, weakly stable cyclops states, and nonstationary regimes in the QIF model.Our reduction approach complements the work by Ratas and Pyragas (2018), which assumed a Lorentzian distribution of the input currents and employed the thermodynamic limit to reduce the QIF model with time-delayed coupling to macroscopic equations that characterize the mean membrane potential, the spiking rate, and the mean synaptic current.The bifurcation analysis of these macroscopic equations, performed in Ratas and Pyragas (2018), revealed alternating parameter regions where the QIF network exhibits macroscopic self-oscillations as a function of input current heterogeneity and time-delayed coupling.While sharing some similarities and goals, our approach is fundamentally different.It reduces finite-size QIF networks with an arbitrary current distribution and arbitrary synaptic activation function G(t/τ) to the finite-size microscopic KS model.This allows for characterizing fine dynamical patterns, including cluster and cyclops states, which might be out of reach for a macroscopic description.In light of this, our reduction approach to an analytically tractable Kuramoto model holds promise in facilitating constructive analysis of rhythmogenesis in QIF networks.By utilizing the reduced KS model, a variety of methods and analytical machinery can be applied to the analysis of collective dynamics in QIF models, including the constructive selection of complex patterns, such as chimeras, whose existence and emergence might be easier to deduce from the reduced Kuramoto model description.Furthermore, the reduction approach holds the potential for extensions to incorporate synaptic adaptation, Hebbian learning, and a complex network structure by employing nodedegree block approximation.

FIGURE 5
FIGURE 5 Dynamical equivalence between the theta neuron (Eqs 3, 4) and KS models (Eq.20), demonstrated via the onset of full synchronization.(A) The evolution of the first |R 1 | (solid curves) and second |R 2 | (dashed curves) order parameters for the theta neuron (green curves) and KS model (red curves), including the transient period.Initial phases θ n , n 1, . . ., N 21 are uniformly distributed over the interval [−π; π].(B) The colors depict the phase differences θ n − θ 21 in the theta neuron model converging to imperfect full synchronization, subject to mismatched parameters η n that are uniformly distributed on the segment [ η − δη/2; η + δη/2], η 2.0, δη 6 × 10 −3 .(C) Comparison between the dynamics of the phase differences θ n − θ 21 for n n 1 17 (thick red curve) and n n 2 20 (thick blue curve) in the theta neuron model and θ n − θ 21 recalculated from phases φ n using the relation (Eq.7) for n n 1 17 (thin cyan curve) and n n 2 20 (thin red curve) in the KS model.Note the perfect alignment of the phase-difference dynamics in the two models.Parameters: q 2, ϰ −0.2π, τ 0.5, ] 20, η 2.0 correspond to point A on Figure 3.